close all;

%% Solver options
tol = 10^-3;
options = odeset('RelTol',tol, 'Stats', 'of');
T = [0:1/49:1];

%% Parameters
paramset = [1, 1, 1;
            1, 0, 100;
            -1, 0, 100;
            -10, 10, 0;
            1, 10000, 0;
            1, 10000, 250;
            10000, 1, 250];
        
func_evals = zeros(size(paramset, 1), 3);
timings = zeros(size(paramset, 1), 3);
         
for experiment = 1:size(paramset, 1)
    params = paramset(experiment, :);
    a = params(1);
    b = params(2);
    c = params(3);
    
    t = 0:1/49:1;
    y_exact = 1 + exp(-b * t) .* sin(c * t) + exp(-a * t);
    [fun, y0, A, F] = stelsel(a, b, c);
    
    %figure();
    eigenvalues = eig(A)
    %plot(eig(A), 'o')
    %title(strcat('Eigenwaarden experminent ', int2str(experiment)));
    
    figure();
    hold on;
    
    colors = ['ro-'; 'go-'; 'yo-'];
    
    i = 0;
    for ode = {@ode45, @ode15s, @ode113}
        i = i + 1;
        time = cputime;
        [x, y, stats] = ode{1}(fun, T, y0, options);
        timings(experiment, i)  = cputime - time;
        
        plot(x, y(:,1), colors(i,:));
        func_evals(experiment, i) = stats(3);
    end
    plot(t, y_exact, 'b-');
    title(strcat('Resultaat experiment ', int2str(experiment)));
end

func_evals
timings